# Script produces analysis data set from source data 
# Last Updated: January 15, 2024

# Load packages 
# tidyverse (Version 1.3.1)
if((!"tidyverse" %in% installed.packages()) | # If package is not installed 
   packageVersion("tidyverse") != "1.3.1"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/tidyverse/tidyverse_1.3.1.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("tidyverse")
} else { # Otherwise, load package
  library("tidyverse")
}

# lubridate (Version 1.8.0)
if((!"lubridate" %in% installed.packages()) | # If package is not installed 
   packageVersion("lubridate") != "1.8.0"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/lubridate/lubridate_1.8.0.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("lubridate")
} else { # Otherwise, load package
  library("lubridate")
}

# gtools (Version 3.9.2.2)
if((!"gtools" %in% installed.packages()) | # If package is not installed 
   packageVersion("gtools") != "3.9.2.2"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/gtools/gtools_3.9.2.2.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("gtools")
} else { # Otherwise, load package
  library("gtools")
}

# Import Raw UNTS Treaty Data 
# Note: This data was scraped from the UNTS website on June 8, 2020 
data <- read.csv("data/treaty_data/unts_source_data.csv") 

# Import Treaty Type (e.g., open multilateral, closed multilateral, etc.)
urls <- read.csv("data/treaty_data/UNTS_urls.csv") 
data <- merge(data, urls, by = "url", all = T) # Merge
rm(urls)

# Clean up UNTS data 
data[data == -999] <- NA # Change missingness code from -999 to NA
data$url <- as.character(data$url) # Coerce vectors to strings
data$Date.of.Notification.Deposit <- as.Date(as.character(data$Date.of.Notification.Deposit), "%d/%m/%Y") # Reformat time variables 
data$Date.of.Effect <- as.Date(as.character(data$Date.of.Effect), "%d/%m/%Y")
data$conclusion_date <- as.Date(as.character(data$conclusion_date), "%d/%m/%Y")
year(data$Date.of.Notification.Deposit)[year(data$Date.of.Notification.Deposit) == 3003] <- 2003  # Fix obvious error

# Compute decision date and year: minimum of date of effect or date of notification 
data$decision_date <- pmin(data$Date.of.Effect, data$Date.of.Notification.Deposit, na.rm = T)
data$decision_date <- as.Date(data$decision_date, origin ="1970-01-01") # Convert back to date
data$decision_year <- as.numeric(format(data$decision_date,'%Y')) # Compute decision year

# Import key for Participants
participant_key <- read.csv("data/treaty_data/participant_key.csv") 
data <- merge(data, participant_key, by = "Participant", all = T) # Merge

# Import key for Actions 
actions_key <- read.csv("data/treaty_data/actions_key.csv") 
actions_key$Action <- as.character(actions_key$Action) # Convert to character
data$Action <- as.character(data$Action)
data <- merge(data, actions_key, by = "Action", all = T) # Merge

# Import key for Issue Areas 
tags_key <- read.csv("data/treaty_data/tags_key.csv") 

# Extract keywords for each issue area 
security <- subset(as.character(tags_key$tags), tags_key$security == 1) # Security
economics <- subset(as.character(tags_key$tags), tags_key$economics == 1) # Economics
human_rights <- subset(as.character(tags_key$tags), tags_key$human_rights == 1) # Human Rights
environment <- subset(as.character(tags_key$tags), tags_key$environment == 1) # Environment

# Create list of each actions subject tags
tags <- as.character(data$subject_terms) # Convert factor to vector
tags[tags == ""] <- 0 # Replace blank values with "0" (necessary b/c missing values in list)
tags <- strsplit(tags, "  ") # Split unique tags within each treaty into unique terms

# Function to identify if each action tag is in a given keyword list
tag_identifier <- function(tag_list, keyword_list){ 
  vals <- lapply(tag_list, function(x) as.numeric(unlist(x) %in% keyword_list)) # Assess if each subject tag is in a list
  vals <- lapply(vals, function(x) x[which.max(x)]) # Extract max value of each list
  vals <- unlist(vals)
  vals
}

# Classify treaties by issue area 
data$economics <- tag_identifier(tags, economics) 
data$security <- tag_identifier(tags, security)
data$human_rights <- tag_identifier(tags, human_rights)
data$environment <- tag_identifier(tags, environment)

# Appendix A.4 ----
tag_freq_fxn <- function(tag_list, keyword_list){
  vals <- lapply(tag_list, function(x) subset(unlist(x), unlist(x) %in% keyword_list)) # Extract economics tags
  vals <- as.data.frame(table(unlist(vals))) # Create frequency table
  vals <- vals[order(vals$Freq, decreasing = T),] # Sort data 
  print(vals)
}

# Frequency of tags at treaty level 
treaty_tags <- data |> group_by(treaty_id) |> summarise(treaty_tags = first(as.character(subject_terms))) # Get tag for each treaty
treaty_tags <- treaty_tags$treaty_tags
treaty_tags[treaty_tags == ""] <- 0 # Replace blank values with "0" (necessary b/c missing values in list)
treaty_tags <- strsplit(treaty_tags, "  ") # Split unique tags within each treaty into unique terms
tag_freq_fxn(treaty_tags, security)
tag_freq_fxn(treaty_tags, economics)
tag_freq_fxn(treaty_tags, environment)
tag_freq_fxn(treaty_tags, human_rights)

# Appendix A Statistics and Table 1 ----
nrow(data) # Total number of actions reported in Appendix A 
summary(data$conclusion_year)
length(unique(data$url)) # Number of unique treaties
length(unique(data$Participant)) # Number of unique participants
length(unique(data$Participant[is.na(data$part_ccode) == F])) # Number of unique ones classified as states
length(unique(data$Action)) # Number of unique actions
sum(data$ratification, na.rm = T) # Number total ratifications
sum(data$withdrawal, na.rm = T) # Number total withdrawals (cooperative and unilateral exits)
nrow(data) - (sum(data$ratification, na.rm = T) + sum(data$withdrawal, na.rm = T)) # Total ambiguous actions
treaty_types <- data |>  # Table 1 in Appendix A
  group_by(url) |> 
  summarize(
    type = first(type),
    action_table = first(action_table)
  )
table(treaty_types$type) # Marginal distributions
table(treaty_types$type[treaty_types$action_table == 1]) # Joint distributions
table(treaty_types$type[treaty_types$action_table == 0]) # Joint distributions
length(unique(data$url)) - sum(table(treaty_types$type)) # Number of treaties missing type classification
rm(treaty_types)

# Also note that data was scraped in June 2020 and there is a delay between when actions 
# are deposited in UNTS and published on its website. Consequently, data quality  
# diminishes beginning in 2019; therefore, I end my study period in 2018.
table(data$conclusion_year[data$conclusion_year > 2015]) # Number treaties concluded each year post-2015
table(data$decision_year[data$decision_year > 2015]) # Number actions each year post-2015

# The UNTS and League of Nations Treaty Series have duplicate registration numbers
# Here we import a list of URLs and prefixes to disambiguate these registration numbers
LoN_II_prefixes <- read.csv("data/treaty_data/LoN_II_treaties_list.csv") # Import list 
data <- merge(data, LoN_II_prefixes, by = "url", all = T) # Merge in LoN data
data$prefix[is.na(data$prefix)] <- "I" # Assign missing values "I" for initial agreement
data <- data |> group_by(registration_number) |> mutate(amendment = as.numeric(conclusion_date != min(conclusion_date))) # Make amendment indicator vairbale
data$prefix[data$amendment == 1] <- "A" # Overwrite "I" with "A" for amendments 
data$prefix <- factor(data$prefix) # Coerce to factor
# Disambiguate LoN and UNTS registration numbers by adding .1 to LoN treaties' reg. numbers
data$registration_number[data$prefix == "LoN"] <- data$registration_number[data$prefix == "LoN"] + 0.1

# Classification of Regimes
regimes <- read.csv("data/treaty_data/regime_ids.csv") # Import treaty-regime key
data <- left_join(data, regimes, by = "treaty_id") # Merge 

# Import COW system data
start_year <- 1945 # Set year to begin analysis/time series
cow_system <- read.csv("data/preprocessing/system2016.csv") # Import COW system data 
cow_system <- subset(cow_system, cow_system$year >= start_year & cow_system$year <= 2018) # subset relevant years
temp_dat <- subset(cow_system, cow_system$year == 2016) # append on 2017 and 2018 years using COW system in 2016
temp_dat$year <- 2017
temp_dat2 <- temp_dat
temp_dat2$year <- 2018
temp_dat <- rbind(temp_dat, temp_dat2)
cow_system <- rbind(cow_system, temp_dat) 
rm(temp_dat, temp_dat2)
cow_system$version <- NULL # Drop version column

# Select years with reliable data  

# Second, note that 1191 treaties are missing conclusion dates
length(unique(data$treaty_id[is.na(data$conclusion_year)]))  
# Note that results are not sensitive to excluding treaties missing conclusion_year values
# For illustration, uncomment the next line and the results will replications
#data <- subset(data, data$conclusion_year <= 2018) 
# However, given considerable missingness, I use the year of the first action taken toward 
# a treaty as its conclusion date in place of these missing values
data <- data |> 
  group_by(treaty_id) |>
  mutate(
    conclusion_year = ifelse(is.na(conclusion_year), min(decision_year), conclusion_year)
  )
length(unique(data$treaty_id[is.na(data$conclusion_year)])) # This leaves 188 treaties missing conclusion dates
data <- subset(data, data$conclusion_year <= 2018) # Select relevant years in treaty data 

# Collapse by treaty_id and compute treaty-level summary statistics 
relevant_treaties <- data |> 
  group_by(treaty_id) |> 
  summarise(
    type = first(agreement_type), # Agreement type 
    url = first(url), # URL
    subject_missing = sum(as.character(subject_terms) == ""), # Missing subject terms
    title = first(title), # Title
    prefix = first(prefix), # Prefix
    min_decision = min(decision_date), # Minimum decision date
    conclusion_year = first(conclusion_year), # Conclusion year
    conclusion_date = first(conclusion_date), # Conclusion date
    registration_number = first(registration_number), # Registration number
    regime_number = first(regime_number), # Regime number 
    core_treaty = max(core_treaty, na.rm = T), # Core treaty
    amendment_code = mean(amendment_code, na.rm = T), # Amendment code 
    n_actions = n(), # Number distinct actions
    n_states = length(unique(part_ccode)), # Number state parties
    n_ratification = sum(ratification), # Number ratifications
    n_withdrawals = sum(withdrawal),  # Number withdrawals
    subject_terms = first(subject_terms), # Subject terms
    environment = max(environment), # Issue areas indicators
    security = max(security),
    human_rights = max(human_rights),
    economics = max(economics)
  ) 

# Zero impute ratifications for treaties without ratifications
relevant_treaties$n_ratification[is.na(relevant_treaties$n_ratification)] <- 0 

# Sort by most significant treaty so that titles are more legible/noteworthy
relevant_treaties <- relevant_treaties[order(-relevant_treaties$n_ratification),] 

# Collapse by regime_number 
multi_treaties <- relevant_treaties |> 
  group_by(regime_number) |> 
  summarize(n = n(), 
            title_of_max_rat = first(title),
            sum_n_ratification = sum(n_ratification), 
            sum_0_ratification = sum(n_ratification == 0), 
            max_n_ratification = max(n_ratification, na.rm = T), 
            n_withdrawals = sum(n_withdrawals, na.rm = T), 
            n_unique_sub = length(unique(subject_terms)),
            environment = sum(environment),
            security = sum(security),
            human_rights = sum(human_rights),
            economics = sum(economics),
            first_conclusion_year = first(conclusion_year),
            min_conclusion_year = min(conclusion_year, na.rm = T),
            median_conclusion_year = median(conclusion_year, na.rm = T),
            mean_conclusion_year = mean(conclusion_year, na.rm = T),
            max_conclusion_year = max(conclusion_year, na.rm = T)
  )

# Select only multilateral treaties (e.g. those with membership > 1)
multi_treaties <- subset(multi_treaties, max_n_ratification > 1) 

# Subset action data 
data <- subset(data, data$regime_number %in% multi_treaties$regime_number) 

# Sort data by year and country so that we loop through the data sequentially, otherwise I might overwrite subsequent actions with earlier actions
sum(is.na(data$decision_date)) # some data are missing decision dates
nrow(subset(data, is.na(data$decision_date) & (data$ratification == 1 | data$withdrawal == 1))) # However only 75 are ratifications or withdrawals 
data <- subset(data, !(is.na(data$decision_date))) # Drop actions missing decision_date
data <- data[order(data$decision_date, data$part_ccode),] # sort data

# Create treaty and regime ratification matrices 
# First, treaty ratification matrices (used to measure indirect treatment and spillover)
data_full <- data # Backup data as data_full
data <- subset(data_full, (data_full$ratification == 1 | data_full$withdrawal == 1)) # Subset to only ratification and withdrawal decisions
data <- subset(data, !(is.na(data$regime_number))) # Drop missing regime numbers

# Set up cooperation matrix
tr <- as.data.frame(matrix(0, nrow = length(cow_system$year), 
                           ncol = (ncol(cow_system) + length(unique(data$treaty_id))))) # Create empty matrix of appropriate dimensions
names <- c(colnames(cow_system), unique(data$treaty_id)) # extract column names from cow_system and treaty id
colnames(tr) <- names # assign column names to matrix
tr[,1:3] <- cow_system # assign cow_system to first four columns of ratification matrix 

# Assigns "ratification" value to tr matrix. 
# Note: When cell == 1, it indicates ratification. 
#       When cell == 0, it indicates withdrawal or non-ratification 
for(i in 1:nrow(data)){
  tr[,as.character(data$treaty_id[i])][tr$year >= data$decision_year[i] & tr$ccode == data$part_ccode[i]] <- data$ratification[i] 
}

# Repeat for withdrawal matrix
wd <- as.data.frame(matrix(0, nrow = length(cow_system$year), 
                           ncol = (ncol(cow_system) + length(unique(data$treaty_id))))) # Create empty matrix of appropriate dimensions
names <- c(colnames(cow_system), unique(data$treaty_id)) # extract column names from cow_system and treaty id
colnames(wd) <- names # assign column names to matrix
wd[,1:3] <- cow_system # assign cow_system to first four columns of ratification matrix. 

# Assigns "withdrawal" value to wd matrix. 
# Note: When cell == 1, it indicates withdrawal having occurred, 
#       When cell == 0, it indicates non-ratification or ongoing ratification  
for(i in 1:nrow(data)){
  wd[,as.character(data$treaty_id[i])][wd$year >= data$decision_year[i] & wd$ccode == data$part_ccode[i]] <- data$withdrawal[i] 
}

# Make tr and wd conformable 
names(wd)[1] <- names(tr)[1] <- "cyear" # rename column
wd$cyear <- tr$cyear <- tr$ccode * 10000 + tr$year # reassign values
idx_names <- names(tr)[1:3] # save index column names to vector
wd <- wd[ , order(names(wd))] # sort columns
tr <- tr[ , order(names(tr))]
tr <- tr |> select(all_of(idx_names), everything()) # move index variables up front
wd <- wd |> select(all_of(idx_names), everything())
sum(names(wd) == names(tr)) == ncol(tr) # check matrices are ordered properly

# Make keys for tr and wd 
tr_key <- tibble(
  treaty_id = names[4:length(names)],
  treaty_names = NA,
  conclusion_year = NA,
  regime_number = NA
)  
for(i in 1:nrow(tr_key)){
  dat_temp <- subset(data, data$treaty_id == tr_key$treaty_id[i])
  tr_key$treaty_names[i] <- dat_temp$title[1]
  tr_key$conclusion_year[i] <-  dat_temp$conclusion_year[1]
  tr_key$regime_number[i] <- dat_temp$regime_number[1]
}

# Save bipartite graphs and key
treaty_ratification_matrices_file_path <- "data/treaty_data/treaty_ratification_matrices.RData"
save(list = c("tr_key","tr", "wd"),  
     file = treaty_ratification_matrices_file_path)

# Second, regime ratification matrices (used to construct panel of regimes)
data <- subset(data, data$core_treaty == 1) # Subset to core_treaty

# Create regime cooperation matrix
tr <- as.data.frame(matrix(0, nrow = length(cow_system$year), 
                           ncol = (ncol(cow_system) + length(unique(data$regime_number))))) # Create empty matrix of appropriate dimensions
names <- c(colnames(cow_system), unique(data$regime_number)) # extract column names from cow_system and treaty id
colnames(tr) <- names # assign column names to matrix
tr[,1:3] <- cow_system # assign cow_system to first four columns of ratification matrix. 

# Assigns "ratification" value to tr matrix. 
# Note: When cell == 1, it indicates ratification. 
#       When cell == 0, it indicates withdrawal or non-ratification 
for(i in 1:nrow(data)){
  tr[,as.character(data$regime_number[i])][tr$year >= data$decision_year[i] & tr$ccode == data$part_ccode[i]] <- data$ratification[i] 
}

# Repeat for regime withdrawal matrix
wd <- as.data.frame(matrix(0, nrow = length(cow_system$year), 
                           ncol = (ncol(cow_system) + length(unique(data$regime_number))))) # Create empty matrix of appropriate dimensions
names <- c(colnames(cow_system), unique(data$regime_number)) # extract column names from cow_system and treaty id
colnames(wd) <- names # assign column names to matrix
wd[,1:3] <- cow_system # assign cow_system to first four columns of ratification matrix. 

# Assigns "withdrawal" value to wd matrix. 
# Note: When cell == 1, it indicates withdrawal having occurred, 
#       When cell == 0, it indicates non-ratification or ongoing ratification 
for(i in 1:nrow(data)){
  wd[,as.character(data$regime_number[i])][wd$year >= data$decision_year[i] & wd$ccode == data$part_ccode[i]] <- data$withdrawal[i] 
}

# Sort tr, wd to make them conformable 
# Make first column cyear
names(wd)[1] <- names(tr)[1] <- "cyear" # rename column
wd$cyear <- tr$cyear <- tr$ccode * 10000 + tr$year # reassign values
idx_names <- names(tr)[1:3] # save index column names to vector
wd <- wd[ , order(names(wd))] # sort columns
tr <- tr[ , order(names(tr))]
tr <- tr |> select(all_of(idx_names), everything()) # move index variables up front
wd <- wd |> select(all_of(idx_names), everything())
sum(names(wd) == names(tr)) == ncol(tr) # check matrices are ordered properly

# Make key regimes
tr_key_regimes <- subset(multi_treaties, multi_treaties$regime_number %in% as.numeric(colnames(tr[,4:ncol(tr)]))) # note one treaty is dropped
tr_regimes <- tr # Disambiguate object names 
wd_regimes <- wd

# Save bipartite graphs and key
regime_ratification_matrices_file_path <- "data/treaty_data/regime_ratification_matrices.RData"
save(list = c("tr_key_regimes","tr_regimes", "wd_regimes"),  
     file = regime_ratification_matrices_file_path)

# Construct Panel of Treaty Regimes
tr_key_regimes <- tr_key_regimes[order(tr_key_regimes$regime_number, tr_key_regimes$min_conclusion_year),] # Sort 
relevant_treaties_multi <- subset(relevant_treaties, # List of all treaties classified as founding treaty or reform
                                  relevant_treaties$type == "Multilateral" & relevant_treaties$amendment_code == 1) 
reform_data <- as.data.frame(as.numeric())

# Loop through regimes, constructing outcome variables and membership indicators
for(i in 1:nrow(tr_key_regimes)){
  # Create panel for each regime during time period and fill in values with relevant data
  dat_temp <- tibble( # Define/initalize key variables for each year in time series
    regime_number = rep(tr_key_regimes$regime_number[i], length(1945:2018)), # Regime number 
    year = 1945:2018, # Year
    title_of_max_rat = tr_key_regimes$title_of_max_rat[i], # Regime name (set to title of most ratified treaty)
    conclusion = as.numeric(tr_key_regimes$min_conclusion_year[i] <= year), # Year of first treaty concluded in regime
    security = NA, human_rights = NA, environment = NA, economics = NA, # Count of treaties in each issue area 
    cumulative_agreements = NA, # Cumulative number of treaties concluded in regime
    n_unique_subject_terms = NA, # Number of unique subject terms
    cumulative_agreements_plus_action_table_amendments = NA, # Cumulative number of treaties concluded in regime plus amendments recorded in action table
    cumulative_unique_subject_terms = NA # Cumulative number of unique subject terms
    )
  relevant_treaties_multi_temp <- subset(relevant_treaties_multi, # Select treaties in given regime
                                         relevant_treaties_multi$regime_number == tr_key_regimes$regime_number[i]) 
  for(y in 1945:2018){ # Loop through each year
    # Get all treaties concluded up until the current year
    concluded_agreements <- subset(relevant_treaties_multi_temp, relevant_treaties_multi_temp$conclusion_year <= y)
    if(nrow(concluded_agreements) > 0){
      # Combine data on all treaties and add to variables for that year
      concluded_agreements_sum <- concluded_agreements |> 
        group_by(regime_number) |> 
        summarize( 
          security = sum(security, na.rm = T), # Sum of treaties within each issue area through year y
          human_rights = sum(human_rights, na.rm = T), 
          environment = sum(environment, na.rm = T), 
          economics = sum(economics, na.rm = T),
          cumulative_agreements = n(), # Total number of treaties through year y
          n_unique_subject_terms = length(unique(subject_terms)) # Total number of unique subject terms through year y
          ) 

      # Check for amendments recorded in action table
      action_amends <- subset(data_full, treaty_id %in% concluded_agreements$treaty_id[i] & 
                                decision_year <= y &
                                Participant == "None" & 
                                Action %in% c("Amendment", "Adoption of Amendments"))
      concluded_agreements_sum$cumulative_agreements_plus_action_table_amendments <- concluded_agreements_sum$cumulative_agreements + nrow(action_amends) 
      concluded_agreements_sum$cumulative_unique_subject_terms <- length(unique(unlist(str_split(unique(concluded_agreements$subject_terms), "  ")))) # Unique number of subject terms
      dat_temp[dat_temp$year == y, 5:12] <- concluded_agreements_sum[,2:9] 
    }
  }
  
  # Compute outcomes
  dat_temp <- dat_temp |> mutate(
    reform_indicator = as.numeric(cumulative_agreements > lag(cumulative_agreements, n = 1)),
    reform_n = as.numeric(cumulative_agreements - lag(cumulative_agreements, n = 1)),
    reform_with_meta_indicator = as.numeric(cumulative_agreements_plus_action_table_amendments > lag(cumulative_agreements_plus_action_table_amendments, n = 1)),
    reform_with_meta_n = as.numeric(cumulative_agreements_plus_action_table_amendments - lag(cumulative_agreements_plus_action_table_amendments, n = 1)),
    new_subject = as.numeric(cumulative_unique_subject_terms != lag(cumulative_unique_subject_terms, n = 1))
    )
  
  # Get ratification data
  tr_temp <- tr_regimes[,c("ccode", "year", as.character(tr_key_regimes$regime_number[i]))] 
  
  # Pivot data
  tr_temp <- tr_temp |> pivot_wider(names_from = ccode, values_from = as.character(tr_key_regimes$regime_number[i]))
  
  # Merge in ratification data 
  dat_temp <- left_join(dat_temp, tr_temp, by = "year")
  
  # Append result to reform_data
  if(nrow(reform_data) == 0){
    reform_data <- dat_temp
  } else {
    reform_data <- rbind(reform_data, dat_temp)
  }
  
}

# Note that new treaties are coded as having new subjects. These should be corrected because subject cannot evolve from nothing. 
reform_data <- reform_data |> group_by(regime_number) |> mutate(new_subject_lag = lag(new_subject, n = 1))
reform_data$new_subject[is.na(reform_data$new_subject_lag)] <- NA

# Import index of withdrawals 
unilateral_exits <- read.csv("data/treaty_data/unilateral_exits.csv")
unilateral_exits <- left_join(unilateral_exits, # Merge in issue area and regime_number from relevant_treaties
                          relevant_treaties[c("treaty_id", "regime_number", 
                                              "economics", "security", "human_rights", "environment")],
                          by = "treaty_id")

# Compute membership size
reform_data$n_members <- rowSums(reform_data[,18:219], na.rm = T)

## National Military Capabilities Data 
nmc <- read.csv("data/preprocessing/NMC-60-abridged.csv") # Load data
# nmc[nmc == - 9] <- NA # Not necessary, because there are no missing CINC values 
nmc_2017 <- subset(nmc, year == 2016) # Use 2016 values for 2017 and 2018 rather than listwise deletion
nmc_2017$year <- 2017
nmc_2018 <- subset(nmc, year == 2016)
nmc_2018$year <- 2018
nmc <- rbind(nmc, nmc_2017, nmc_2018)
nmc$cyear <- nmc$ccode*10000 + nmc$year # Generate keys for merging
unilateral_exits$cyear <- unilateral_exits$part_ccode*10000 + unilateral_exits$decision_year 
unilateral_exits <- left_join(unilateral_exits, nmc[c("cyear", "cinc")], by = "cyear") # Merge

# Relative Power of Exiter w.r.t. Treaty Members
unilateral_exits$relative_power_exiter <- NA
load(file = treaty_ratification_matrices_file_path) # Import treaty data
for(i in 1:nrow(unilateral_exits)){ # Loop through withdrawals
  if(unilateral_exits$treaty_id[i] %in% colnames(tr)){
    members_temp <- subset(tr$cyear, # Get cyear for treaty members 
                           tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                             tr$year ==  unilateral_exits$decision_year[i])
    member_cincs <- nmc$cinc[nmc$cyear %in% members_temp] # Get treaty member cinc  
    exiter_cinc <- nmc$cinc[nmc$cyear == unilateral_exits$cyear[i]] # Get exiter cinc
    mu_members <- mean(member_cincs, na.rm = T) # Get treaty member mean
    power_sum <- sum(c(exiter_cinc, member_cincs), na.rm = T) # Get treaty total power
    if(length(exiter_cinc) > 0 & length(member_cincs) > 0){ # If cinc scores are available
      unilateral_exits$relative_power_exiter[i] <- exiter_cinc/power_sum # Get exiter's proportion of total power 
    }
  }
}

power_cutoff_quantile <- 0.5 # Use median as cutoff value for main analysis
unilateral_exits$rel_power <- as.numeric(unilateral_exits$relative_power_exiter > quantile(unilateral_exits$relative_power_exiter, probs = power_cutoff_quantile, na.rm = T))
unilateral_exits$rel_power_tertile <- as.numeric(quantcut(unilateral_exits$relative_power_exiter, 3))

# Relative Power of Exiter w.r.t. All Exiters
unilateral_exits$power <- as.numeric(unilateral_exits$cinc > quantile(unilateral_exits$cinc, probs = power_cutoff_quantile, na.rm = T))
unilateral_exits$power_tertile <- as.numeric(quantcut(unilateral_exits$cinc, 3))

# If two different states exit a treaty in a regime in the same year, then
# use the greater of the two states power to define treatment
unilateral_exits <- unilateral_exits |> 
  group_by(decision_year, regime_number) |>
  mutate(
    power = max(power, na.rm = T),
    rel_power = max(rel_power, na.rm = T),
    power_tertile = max(power_tertile, na.rm = T),
    rel_power_tertile = max(rel_power_tertile, na.rm = T)   
  )

# Correct -Inf created by missing regime_numbers (though this doesn't matter)
unilateral_exits$rel_power[unilateral_exits$rel_power < 0] <- NA
unilateral_exits$rel_power_tertile[unilateral_exits$rel_power_tertile < 0] <- NA

# Create treatment variable in reform_data
reform_data$direct_exit_initial <- 0 # Treatment variable denoting that this regime had a state exit unilaterally from a component treaty in year t
reform_data$indirect_exit_initial <- 0 # Treatment variable denoting that a state in this regime exited unilaterally from any component treaty in year t
reform_data$total_affected_members <- 0 # Treatment variable denoting that a state in this regime exited unilaterally from any component treaty in year t and there is at least one other member of the affected treaty in this regime

reform_data$direct_exit_initial_power <- 0 
reform_data$indirect_exit_initial_power <- 0 
reform_data$total_affected_members_power <- 0 

reform_data$direct_exit_initial_not_power <- 0 
reform_data$indirect_exit_initial_not_power <- 0 
reform_data$total_affected_members_not_power <- 0 

reform_data$direct_exit_initial_rel_power <- 0 
reform_data$indirect_exit_initial_rel_power <- 0 
reform_data$total_affected_members_rel_power <- 0 

reform_data$direct_exit_initial_not_rel_power <- 0 
reform_data$indirect_exit_initial_not_rel_power <- 0 
reform_data$total_affected_members_not_rel_power <- 0 

reform_data$direct_exit_initial_power_tertile1 <- 0 
reform_data$indirect_exit_initial_power_tertile1 <- 0 
reform_data$total_affected_members_power_tertile1 <- 0 

reform_data$direct_exit_initial_power_tertile2 <- 0 
reform_data$indirect_exit_initial_power_tertile2 <- 0 
reform_data$total_affected_members_power_tertile2 <- 0 

reform_data$direct_exit_initial_power_tertile3 <- 0 
reform_data$indirect_exit_initial_power_tertile3 <- 0 
reform_data$total_affected_members_power_tertile3 <- 0 

reform_data$direct_exit_initial_rel_power_tertile1 <- 0 
reform_data$indirect_exit_initial_rel_power_tertile1 <- 0 
reform_data$total_affected_members_rel_power_tertile1 <- 0 

reform_data$direct_exit_initial_rel_power_tertile2 <- 0 
reform_data$indirect_exit_initial_rel_power_tertile2 <- 0 
reform_data$total_affected_members_rel_power_tertile2 <- 0 

reform_data$direct_exit_initial_rel_power_tertile3 <- 0 
reform_data$indirect_exit_initial_rel_power_tertile3 <- 0 
reform_data$total_affected_members_rel_power_tertile3 <- 0 

# Loop through withdrawals
for(i in 1:nrow(unilateral_exits)){

  #' The following sections of code compute indicator variables for regimes that are 
  #' directly and indirectly treated and the number of states in each regime 
  #' that were party to the regime where withdrawal occurred. This last variable, 
  #' "total_affected_members", is used to compute overlap below. Each section
  #' has the same code except for the initial "if" statement that determines
  #' which withdrawal from unilateral_exits to define each treatment. 

  # 1. All initial unilateral exits
  if(is.na(unilateral_exits$regime_number[i]) == F){ # Define direct and indirect treatment indicator variables 
    reform_data$direct_exit_initial[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1 # Equals 1 if withdrawal happened in regime's constituent treaty in that year
    reform_data$indirect_exit_initial[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1 # Equals 1 if treaty member withdrew from another agreement in that year (below I fix cases where there is are no other overlapping members)
  }
  if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){   # Compute number of members in each regime that were party to the treaty where withdrawal occurred
    members_temp <- subset(tr$ccode, # Make list of other state parties
                           tr[, as.character(unilateral_exits$treaty_id[i])] == 1 & 
                             tr$year == unilateral_exits$decision_year[i])
    members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) # Exclude exiting state
    reform_data$total_affected_members[reform_data$year== unilateral_exits$decision_year[i]] <- # How many states in this regime were in the treaty where withdrawal occurred?
      rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = F) + 
      reform_data$total_affected_members[reform_data$year== unilateral_exits$decision_year[i]] # Sum must be cumulative (e.g., sum over previous withdrawals in loop) 
  }
  
  # 2. Absolute Power (binary)
  if(unilateral_exits$power[i] == 1 & is.na(unilateral_exits$power[i]) == F){ 
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_power[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_power[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_power[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_power[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 3. Not Absolute Power (binary)
  if(unilateral_exits$power[i] == 0 & is.na(unilateral_exits$power[i]) == F){ 
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_not_power[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_not_power[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_not_power[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_not_power[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 4. Absolute Power (Tertile 1)
  if(unilateral_exits$power_tertile[i] == 1 & is.na(unilateral_exits$power_tertile[i]) == F){
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_power_tertile1[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_power_tertile1[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode,
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_power_tertile1[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_power_tertile1[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 5. Absolute Power (Tertile 2)
  if(unilateral_exits$power_tertile[i] == 2 & is.na(unilateral_exits$power_tertile[i]) == F){
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_power_tertile2[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_power_tertile2[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_power_tertile2[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_power_tertile2[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 6. Absolute Power (Tertile 3)
  if(unilateral_exits$power_tertile[i] == 3 & is.na(unilateral_exits$power_tertile[i]) == F){
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_power_tertile3[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_power_tertile3[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_power_tertile3[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_power_tertile3[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 7. Relative Power (binary)
  if(unilateral_exits$rel_power[i] == 1 & is.na(unilateral_exits$rel_power[i]) == F){ 
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_rel_power[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_rel_power[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_rel_power[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_rel_power[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 8. Not Relative Power (binary)
  if(unilateral_exits$rel_power[i] == 0 & is.na(unilateral_exits$rel_power[i]) == F){ 
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_not_rel_power[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_not_rel_power[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_not_rel_power[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_not_rel_power[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 9. Relative Power (Tertile 1)
  if(unilateral_exits$rel_power_tertile[i] == 1 & is.na(unilateral_exits$rel_power_tertile[i]) == F){
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_rel_power_tertile1[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_rel_power_tertile1[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_rel_power_tertile1[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_rel_power_tertile1[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 10. Relative Power (Tertile 2)
  if(unilateral_exits$rel_power_tertile[i] == 2 & is.na(unilateral_exits$rel_power_tertile[i]) == F){
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_rel_power_tertile2[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_rel_power_tertile2[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_rel_power_tertile2[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_rel_power_tertile2[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
  
  # 11. Relative Power (Tertile 3)
  if(unilateral_exits$rel_power_tertile[i] == 3 & is.na(unilateral_exits$rel_power_tertile[i]) == F){
    if(is.na(unilateral_exits$regime_number[i]) == F){
      reform_data$direct_exit_initial_rel_power_tertile3[reform_data$regime_number == unilateral_exits$regime_number[i] & reform_data$year == unilateral_exits$decision_year[i]] <- 1
      reform_data$indirect_exit_initial_rel_power_tertile3[reform_data[,as.character(unilateral_exits$part_ccode[i])] == 1 & reform_data$year == unilateral_exits$decision_year[i]] <- 1
    }
    if(as.character(unilateral_exits$treaty_id[i]) %in% colnames(tr)){
      members_temp <- subset(tr$ccode, 
                             tr[,as.character(unilateral_exits$treaty_id[i])] == 1 & 
                               tr$year ==  unilateral_exits$decision_year[i])
      members_temp <- subset(members_temp, !(members_temp %in% unilateral_exits$part_ccode[i])) 
      reform_data$total_affected_members_rel_power_tertile3[reform_data$year== unilateral_exits$decision_year[i]] <- 
        rowSums(subset(reform_data[,as.character(members_temp)], reform_data$year== unilateral_exits$decision_year[i]), na.rm = T) + 
        reform_data$total_affected_members_rel_power_tertile3[reform_data$year== unilateral_exits$decision_year[i]] 
    }
  }
}

# Measure of average overlap
# First, Compute how many withdrawals per year
n_wd_per_year <- as.data.frame(table(unilateral_exits$decision_year))
colnames(n_wd_per_year) <- c("year", "n_wd")
n_wd_per_year <- unilateral_exits |> 
  group_by(decision_year) |> summarize(
    n_wd = n(),

    n_wd_power = sum(power[is.na(power) == F] == 1),
    n_wd_not_power = sum(power[is.na(power) == F] == 0),
    
    n_wd_power_tertile1 = sum(power_tertile[is.na(power_tertile) == F] == 1),
    n_wd_power_tertile2 = sum(power_tertile[is.na(power_tertile) == F] == 2),
    n_wd_power_tertile3 = sum(power_tertile[is.na(power_tertile) == F] == 3),
    
    n_wd_rel_power = sum(rel_power[is.na(rel_power) == F] == 1),
    n_wd_not_rel_power = sum(rel_power[is.na(rel_power) == F] == 0),
    
    n_wd_rel_power_tertile1 = sum(rel_power_tertile[is.na(rel_power_tertile) == F] == 1),
    n_wd_rel_power_tertile2 = sum(rel_power_tertile[is.na(rel_power_tertile) == F] == 2),
    n_wd_rel_power_tertile3 = sum(rel_power_tertile[is.na(rel_power_tertile) == F] == 3),
  )

# Second, multiply n_members by the number of withdrawals per year
reform_data$n_members_wd_weighted <- 0 # Initialize variables
reform_data$n_members_wd_weighted_power <- 0
reform_data$n_members_wd_weighted_not_power <- 0
reform_data$n_members_wd_weighted_power_tertile1  <- 0
reform_data$n_members_wd_weighted_power_tertile2  <- 0
reform_data$n_members_wd_weighted_power_tertile3  <- 0
reform_data$n_members_wd_weighted_rel_power <- 0
reform_data$n_members_wd_weighted_not_rel_power <- 0
reform_data$n_members_wd_weighted_rel_power_tertile1  <- 0
reform_data$n_members_wd_weighted_rel_power_tertile2  <- 0
reform_data$n_members_wd_weighted_rel_power_tertile3  <- 0
for(i in 1:nrow(n_wd_per_year)){ # Multiply
  reform_data$n_members_wd_weighted[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd[i]
  
  reform_data$n_members_wd_weighted_power[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_power[i]
  
  reform_data$n_members_wd_weighted_not_power[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_not_power[i]
  
  reform_data$n_members_wd_weighted_power_tertile1[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_power_tertile1[i]
  
  reform_data$n_members_wd_weighted_power_tertile2[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_power_tertile2[i]

  reform_data$n_members_wd_weighted_power_tertile3[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_power_tertile3[i]
  
  reform_data$n_members_wd_weighted_rel_power[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_rel_power[i]
  
  reform_data$n_members_wd_weighted_not_rel_power[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_not_rel_power[i]
  
  reform_data$n_members_wd_weighted_rel_power_tertile1[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_rel_power_tertile1[i]
  
  reform_data$n_members_wd_weighted_rel_power_tertile2[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_rel_power_tertile2[i]
  
  reform_data$n_members_wd_weighted_rel_power_tertile3[reform_data$year ==  n_wd_per_year$decision_year[i]] <- 
    reform_data$n_members[reform_data$year ==  n_wd_per_year$decision_year[i]]  * n_wd_per_year$n_wd_rel_power_tertile3[i]
}

# Third, compute average overlap 
# = total number of members exposed to treatment in year /(number members  in year * number of withdrawals in year)
reform_data <- reform_data |> 
  mutate(overlap = total_affected_members/(n_members_wd_weighted), 

         overlap_power = total_affected_members_power/(n_members_wd_weighted_power),
         overlap_not_power = total_affected_members_not_power/(n_members_wd_weighted_not_power),
         
         overlap_power_tertile1 = total_affected_members_power_tertile1/(n_members_wd_weighted_power_tertile1),
         overlap_power_tertile2 = total_affected_members_power_tertile2/(n_members_wd_weighted_power_tertile2),
         overlap_power_tertile3 = total_affected_members_power_tertile3/(n_members_wd_weighted_power_tertile3),
         
         overlap_rel_power = total_affected_members_rel_power/(n_members_wd_weighted_rel_power),
         overlap_not_rel_power = total_affected_members_not_rel_power/(n_members_wd_weighted_not_rel_power),

         overlap_rel_power_tertile1 = total_affected_members_rel_power_tertile1/(n_members_wd_weighted_rel_power_tertile1),
         overlap_rel_power_tertile2 = total_affected_members_rel_power_tertile2/(n_members_wd_weighted_rel_power_tertile2),
         overlap_rel_power_tertile3 = total_affected_members_rel_power_tertile3/(n_members_wd_weighted_rel_power_tertile3),
  ) 

# Two issues to correct:
# First, select only regimes with more than one state party
reform_data <-  subset(reform_data, n_members > 1)  

# Second, fix NAs created by dividing 0/0
table(reform_data$total_affected_members[is.na(reform_data$overlap)]) # Note all NAs have 0 affected members
reform_data$overlap[is.na(reform_data$overlap)] <- 0 # Replace with 0 
reform_data$overlap_power[is.na(reform_data$overlap_power)] <- 0  
reform_data$overlap_not_power[is.na(reform_data$overlap_not_power)] <- 0 
reform_data$overlap_power_tertile1[is.na(reform_data$overlap_power_tertile1)] <- 0  
reform_data$overlap_power_tertile2[is.na(reform_data$overlap_power_tertile2)] <- 0  
reform_data$overlap_power_tertile3[is.na(reform_data$overlap_power_tertile3)] <- 0  
reform_data$overlap_rel_power[is.na(reform_data$overlap_rel_power)] <- 0  
reform_data$overlap_not_rel_power[is.na(reform_data$overlap_not_rel_power)] <- 0  
reform_data$overlap_rel_power_tertile1[is.na(reform_data$overlap_rel_power_tertile1)] <- 0  
reform_data$overlap_rel_power_tertile2[is.na(reform_data$overlap_rel_power_tertile2)] <- 0  
reform_data$overlap_rel_power_tertile3[is.na(reform_data$overlap_rel_power_tertile3)] <- 0  

# Compute years since last reform
reform_data$reform_indicator0 <- reform_data$reform_indicator # Define new outcome variable with NAs zero imputed
reform_data$reform_indicator0[is.na(reform_data$reform_indicator0)] <- 0 # Zero impute NAs 
reform_data <- reform_data |> 
  group_by(regime_number) |> 
  arrange(year) |> 
  mutate(
    last_reform = (ave(reform_indicator0, cumsum(c(F, diff(reform_indicator0) < 0)), FUN=seq_along))
  )

# Coerce regime_number to integer
reform_data$regime_number_int <- as.numeric(factor(reform_data$regime_number))

# Address -5 in reform_with_meta_n (issue is only with one observation)
reform_data$reform_with_meta_n[reform_data$reform_with_meta_n < 0] <- 0

# Regime must have at least 1 overlapping member aside from exiter to qualify as indirectly treated
reform_data$indirect_exit_initial[reform_data$overlap <= 0] <- 0 
reform_data$indirect_exit_initial_power[reform_data$overlap_power <= 0] <- 0 
reform_data$indirect_exit_initial_not_power[reform_data$overlap_not_power <= 0] <- 0 
reform_data$indirect_exit_initial_power_tertile1[reform_data$overlap_power_tertile1 <= 0] <- 0 
reform_data$indirect_exit_initial_power_tertile2[reform_data$overlap_power_tertile2 <= 0] <- 0 
reform_data$indirect_exit_initial_power_tertile3[reform_data$overlap_power_tertile3 <= 0] <- 0 
reform_data$indirect_exit_initial_rel_power_tertile1[reform_data$overlap_rel_power_tertile1 <= 0] <- 0 
reform_data$indirect_exit_initial_rel_power_tertile2[reform_data$overlap_rel_power_tertile2 <= 0] <- 0 
reform_data$indirect_exit_initial_rel_power_tertile3[reform_data$overlap_rel_power_tertile3 <= 0] <- 0 

# Lead variables for placebo analyses 
reform_data <- reform_data |> 
  group_by(regime_number) |> 
  mutate(
    
    direct_exit_initial_placebo_1 = lead(direct_exit_initial, n = 1, order_by = year),
    direct_exit_initial_placebo_2 = lead(direct_exit_initial, n = 2, order_by = year),
    direct_exit_initial_placebo_3 = lead(direct_exit_initial, n = 3, order_by = year),
    direct_exit_initial_placebo_4 = lead(direct_exit_initial, n = 4, order_by = year),
    direct_exit_initial_placebo_5 = lead(direct_exit_initial, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_1 = lead(indirect_exit_initial, n = 1, order_by = year),
    indirect_exit_initial_placebo_2 = lead(indirect_exit_initial, n = 2, order_by = year),
    indirect_exit_initial_placebo_3 = lead(indirect_exit_initial, n = 3, order_by = year),
    indirect_exit_initial_placebo_4 = lead(indirect_exit_initial, n = 4, order_by = year),
    indirect_exit_initial_placebo_5 = lead(indirect_exit_initial, n = 5, order_by = year),
    
    direct_exit_initial_placebo_power_1 = lead(direct_exit_initial_power, n = 1, order_by = year),
    direct_exit_initial_placebo_power_2 = lead(direct_exit_initial_power, n = 2, order_by = year),
    direct_exit_initial_placebo_power_3 = lead(direct_exit_initial_power, n = 3, order_by = year),
    direct_exit_initial_placebo_power_4 = lead(direct_exit_initial_power, n = 4, order_by = year),
    direct_exit_initial_placebo_power_5 = lead(direct_exit_initial_power, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_power_1 = lead(indirect_exit_initial_power, n = 1, order_by = year),
    indirect_exit_initial_placebo_power_2 = lead(indirect_exit_initial_power, n = 2, order_by = year),
    indirect_exit_initial_placebo_power_3 = lead(indirect_exit_initial_power, n = 3, order_by = year),
    indirect_exit_initial_placebo_power_4 = lead(indirect_exit_initial_power, n = 4, order_by = year),
    indirect_exit_initial_placebo_power_5 = lead(indirect_exit_initial_power, n = 5, order_by = year),
    
    direct_exit_initial_placebo_not_power_1 = lead(direct_exit_initial_not_power, n = 1, order_by = year),
    direct_exit_initial_placebo_not_power_2 = lead(direct_exit_initial_not_power, n = 2, order_by = year),
    direct_exit_initial_placebo_not_power_3 = lead(direct_exit_initial_not_power, n = 3, order_by = year),
    direct_exit_initial_placebo_not_power_4 = lead(direct_exit_initial_not_power, n = 4, order_by = year),
    direct_exit_initial_placebo_not_power_5 = lead(direct_exit_initial_not_power, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_not_power_1 = lead(indirect_exit_initial_not_power, n = 1, order_by = year),
    indirect_exit_initial_placebo_not_power_2 = lead(indirect_exit_initial_not_power, n = 2, order_by = year),
    indirect_exit_initial_placebo_not_power_3 = lead(indirect_exit_initial_not_power, n = 3, order_by = year),
    indirect_exit_initial_placebo_not_power_4 = lead(indirect_exit_initial_not_power, n = 4, order_by = year),
    indirect_exit_initial_placebo_not_power_5 = lead(indirect_exit_initial_not_power, n = 5, order_by = year),
    
    direct_exit_initial_placebo_power_tertile1_1 = lead(direct_exit_initial_power_tertile1, n = 1, order_by = year),
    direct_exit_initial_placebo_power_tertile1_2 = lead(direct_exit_initial_power_tertile1, n = 2, order_by = year),
    direct_exit_initial_placebo_power_tertile1_3 = lead(direct_exit_initial_power_tertile1, n = 3, order_by = year),
    direct_exit_initial_placebo_power_tertile1_4 = lead(direct_exit_initial_power_tertile1, n = 4, order_by = year),
    direct_exit_initial_placebo_power_tertile1_5 = lead(direct_exit_initial_power_tertile1, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_power_tertile1_1 = lead(indirect_exit_initial_power_tertile1, n = 1, order_by = year),
    indirect_exit_initial_placebo_power_tertile1_2 = lead(indirect_exit_initial_power_tertile1, n = 2, order_by = year),
    indirect_exit_initial_placebo_power_tertile1_3 = lead(indirect_exit_initial_power_tertile1, n = 3, order_by = year),
    indirect_exit_initial_placebo_power_tertile1_4 = lead(indirect_exit_initial_power_tertile1, n = 4, order_by = year),
    indirect_exit_initial_placebo_power_tertile1_5 = lead(indirect_exit_initial_power_tertile1, n = 5, order_by = year),
    
    direct_exit_initial_placebo_power_tertile2_1 = lead(direct_exit_initial_power_tertile2, n = 1, order_by = year),
    direct_exit_initial_placebo_power_tertile2_2 = lead(direct_exit_initial_power_tertile2, n = 2, order_by = year),
    direct_exit_initial_placebo_power_tertile2_3 = lead(direct_exit_initial_power_tertile2, n = 3, order_by = year),
    direct_exit_initial_placebo_power_tertile2_4 = lead(direct_exit_initial_power_tertile2, n = 4, order_by = year),
    direct_exit_initial_placebo_power_tertile2_5 = lead(direct_exit_initial_power_tertile2, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_power_tertile2_1 = lead(indirect_exit_initial_power_tertile2, n = 1, order_by = year),
    indirect_exit_initial_placebo_power_tertile2_2 = lead(indirect_exit_initial_power_tertile2, n = 2, order_by = year),
    indirect_exit_initial_placebo_power_tertile2_3 = lead(indirect_exit_initial_power_tertile2, n = 3, order_by = year),
    indirect_exit_initial_placebo_power_tertile2_4 = lead(indirect_exit_initial_power_tertile2, n = 4, order_by = year),
    indirect_exit_initial_placebo_power_tertile2_5 = lead(indirect_exit_initial_power_tertile2, n = 5, order_by = year),
    
    direct_exit_initial_placebo_power_tertile3_1 = lead(direct_exit_initial_power_tertile3, n = 1, order_by = year),
    direct_exit_initial_placebo_power_tertile3_2 = lead(direct_exit_initial_power_tertile3, n = 2, order_by = year),
    direct_exit_initial_placebo_power_tertile3_3 = lead(direct_exit_initial_power_tertile3, n = 3, order_by = year),
    direct_exit_initial_placebo_power_tertile3_4 = lead(direct_exit_initial_power_tertile3, n = 4, order_by = year),
    direct_exit_initial_placebo_power_tertile3_5 = lead(direct_exit_initial_power_tertile3, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_power_tertile3_1 = lead(indirect_exit_initial_power_tertile3, n = 1, order_by = year),
    indirect_exit_initial_placebo_power_tertile3_2 = lead(indirect_exit_initial_power_tertile3, n = 2, order_by = year),
    indirect_exit_initial_placebo_power_tertile3_3 = lead(indirect_exit_initial_power_tertile3, n = 3, order_by = year),
    indirect_exit_initial_placebo_power_tertile3_4 = lead(indirect_exit_initial_power_tertile3, n = 4, order_by = year),
    indirect_exit_initial_placebo_power_tertile3_5 = lead(indirect_exit_initial_power_tertile3, n = 5, order_by = year),
    
    direct_exit_initial_placebo_rel_power_1 = lead(direct_exit_initial_rel_power, n = 1, order_by = year),
    direct_exit_initial_placebo_rel_power_2 = lead(direct_exit_initial_rel_power, n = 2, order_by = year),
    direct_exit_initial_placebo_rel_power_3 = lead(direct_exit_initial_rel_power, n = 3, order_by = year),
    direct_exit_initial_placebo_rel_power_4 = lead(direct_exit_initial_rel_power, n = 4, order_by = year),
    direct_exit_initial_placebo_rel_power_5 = lead(direct_exit_initial_rel_power, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_rel_power_1 = lead(indirect_exit_initial_rel_power, n = 1, order_by = year),
    indirect_exit_initial_placebo_rel_power_2 = lead(indirect_exit_initial_rel_power, n = 2, order_by = year),
    indirect_exit_initial_placebo_rel_power_3 = lead(indirect_exit_initial_rel_power, n = 3, order_by = year),
    indirect_exit_initial_placebo_rel_power_4 = lead(indirect_exit_initial_rel_power, n = 4, order_by = year),
    indirect_exit_initial_placebo_rel_power_5 = lead(indirect_exit_initial_rel_power, n = 5, order_by = year),
    
    direct_exit_initial_placebo_not_rel_power_1 = lead(direct_exit_initial_not_rel_power, n = 1, order_by = year),
    direct_exit_initial_placebo_not_rel_power_2 = lead(direct_exit_initial_not_rel_power, n = 2, order_by = year),
    direct_exit_initial_placebo_not_rel_power_3 = lead(direct_exit_initial_not_rel_power, n = 3, order_by = year),
    direct_exit_initial_placebo_not_rel_power_4 = lead(direct_exit_initial_not_rel_power, n = 4, order_by = year),
    direct_exit_initial_placebo_not_rel_power_5 = lead(direct_exit_initial_not_rel_power, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_not_rel_power_1 = lead(indirect_exit_initial_not_rel_power, n = 1, order_by = year),
    indirect_exit_initial_placebo_not_rel_power_2 = lead(indirect_exit_initial_not_rel_power, n = 2, order_by = year),
    indirect_exit_initial_placebo_not_rel_power_3 = lead(indirect_exit_initial_not_rel_power, n = 3, order_by = year),
    indirect_exit_initial_placebo_not_rel_power_4 = lead(indirect_exit_initial_not_rel_power, n = 4, order_by = year),
    indirect_exit_initial_placebo_not_rel_power_5 = lead(indirect_exit_initial_not_rel_power, n = 5, order_by = year),
    
    direct_exit_initial_placebo_rel_power_tertile1_1 = lead(direct_exit_initial_rel_power_tertile1, n = 1, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile1_2 = lead(direct_exit_initial_rel_power_tertile1, n = 2, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile1_3 = lead(direct_exit_initial_rel_power_tertile1, n = 3, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile1_4 = lead(direct_exit_initial_rel_power_tertile1, n = 4, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile1_5 = lead(direct_exit_initial_rel_power_tertile1, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_rel_power_tertile1_1 = lead(indirect_exit_initial_rel_power_tertile1, n = 1, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile1_2 = lead(indirect_exit_initial_rel_power_tertile1, n = 2, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile1_3 = lead(indirect_exit_initial_rel_power_tertile1, n = 3, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile1_4 = lead(indirect_exit_initial_rel_power_tertile1, n = 4, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile1_5 = lead(indirect_exit_initial_rel_power_tertile1, n = 5, order_by = year),
    
    direct_exit_initial_placebo_rel_power_tertile2_1 = lead(direct_exit_initial_rel_power_tertile2, n = 1, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile2_2 = lead(direct_exit_initial_rel_power_tertile2, n = 2, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile2_3 = lead(direct_exit_initial_rel_power_tertile2, n = 3, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile2_4 = lead(direct_exit_initial_rel_power_tertile2, n = 4, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile2_5 = lead(direct_exit_initial_rel_power_tertile2, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_rel_power_tertile2_1 = lead(indirect_exit_initial_rel_power_tertile2, n = 1, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile2_2 = lead(indirect_exit_initial_rel_power_tertile2, n = 2, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile2_3 = lead(indirect_exit_initial_rel_power_tertile2, n = 3, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile2_4 = lead(indirect_exit_initial_rel_power_tertile2, n = 4, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile2_5 = lead(indirect_exit_initial_rel_power_tertile2, n = 5, order_by = year),
    
    direct_exit_initial_placebo_rel_power_tertile3_1 = lead(direct_exit_initial_rel_power_tertile3, n = 1, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile3_2 = lead(direct_exit_initial_rel_power_tertile3, n = 2, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile3_3 = lead(direct_exit_initial_rel_power_tertile3, n = 3, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile3_4 = lead(direct_exit_initial_rel_power_tertile3, n = 4, order_by = year),
    direct_exit_initial_placebo_rel_power_tertile3_5 = lead(direct_exit_initial_rel_power_tertile3, n = 5, order_by = year),
    
    indirect_exit_initial_placebo_rel_power_tertile3_1 = lead(indirect_exit_initial_rel_power_tertile3, n = 1, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile3_2 = lead(indirect_exit_initial_rel_power_tertile3, n = 2, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile3_3 = lead(indirect_exit_initial_rel_power_tertile3, n = 3, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile3_4 = lead(indirect_exit_initial_rel_power_tertile3, n = 4, order_by = year),
    indirect_exit_initial_placebo_rel_power_tertile3_5 = lead(indirect_exit_initial_rel_power_tertile3, n = 5, order_by = year),
  )

# Import regime ideal point values
load("data/preprocessing/AgreementScoresAll_Jul2023.Rdata") # Import years of UNGA sessions
dfAgree <- dfAgree |> group_by(session.x) |> summarize(year = mean(year, na.rm = T))
colnames(dfAgree) <- c("session", "year")
ideal_points <- read.csv("data/preprocessing/IdealpointestimatesAll_Jul2023.csv") # Import UNGA Ideal Point data 
ideal_points <- left_join(ideal_points, dfAgree, by = "session") # Merge

# Initialize variables
reform_data$pref_mean <- NA
reform_data$pref_sd <- NA

for(i in 1:nrow(reform_data)){
  members_temp <- colnames(reform_data)[which(reform_data[i, 18:219] == 1) + 17] # Get treaty members
  year_temp <- reform_data$year[i] # Get year
  ideal_points_temp <- subset(ideal_points$IdealPointAll, # Extract vector of ideal points
                              ideal_points$year == year_temp & 
                              ideal_points$ccode %in% members_temp)
  if(length(ideal_points_temp) > 1){ # If more than one member, compute values
    reform_data$pref_mean[i] <- mean(ideal_points_temp, na.rm = T)
    reform_data$pref_sd[i] <- sd(ideal_points_temp, na.rm = T)
  }
  if(i %% 1000 == 0){print(i)} # Print every 1000th i to check progress
}

# Log reform_n and n_members
reform_data <- reform_data |> 
  mutate(
    reform_n_log = log(reform_n + 1),
    n_members_log = log(n_members + 1)
    )

# Coerce year to integer for PanelMatch 
reform_data$year <- as.integer(reform_data$year)

# Select variables used in analysis
analysis_dat_vars <- c(
  
  # Unit Index Variables
  "title_of_max_rat",
  "year",
  "regime_number",
  "regime_number_int",

  # Outcome
  "reform_n",
  
  # Covariates
  "n_members_log",  
  "last_reform", 
  "cumulative_agreements",
  "pref_mean", 
  "pref_sd",
  "human_rights", 
  "environment",
  "security", 
  "economics",
  
  # Treatment Variables
  ## Direct Effect
  "direct_exit_initial", 
  "direct_exit_initial_placebo_1",  
  "direct_exit_initial_placebo_2", 
  "direct_exit_initial_placebo_3", 
  "direct_exit_initial_placebo_4",  
  "direct_exit_initial_placebo_5", 
  
  ## Direct Effect (Relative Power == 1) 
  "direct_exit_initial_rel_power", 
  "direct_exit_initial_placebo_rel_power_1",  
  "direct_exit_initial_placebo_rel_power_2", 
  "direct_exit_initial_placebo_rel_power_3",  
  "direct_exit_initial_placebo_rel_power_4", 
  "direct_exit_initial_placebo_rel_power_5", 
  
  ## Direct Effect (Relative Power == 0) 
  "direct_exit_initial_not_rel_power", 
  "direct_exit_initial_placebo_not_rel_power_1",  
  "direct_exit_initial_placebo_not_rel_power_2", 
  "direct_exit_initial_placebo_not_rel_power_3",  
  "direct_exit_initial_placebo_not_rel_power_4", 
  "direct_exit_initial_placebo_not_rel_power_5", 
  
  ## Direct Effect (Power Tertile 1) 
  "direct_exit_initial_rel_power_tertile1", 
  "direct_exit_initial_placebo_rel_power_tertile1_1", 
  "direct_exit_initial_placebo_rel_power_tertile1_2", 
  "direct_exit_initial_placebo_rel_power_tertile1_3", 
  "direct_exit_initial_placebo_rel_power_tertile1_4",
  "direct_exit_initial_placebo_rel_power_tertile1_5", 
  
  ## Direct Effect (Power Tertile 2) 
  "direct_exit_initial_rel_power_tertile2", 
  "direct_exit_initial_placebo_rel_power_tertile2_1",  
  "direct_exit_initial_placebo_rel_power_tertile2_2",  
  "direct_exit_initial_placebo_rel_power_tertile2_3",  
  "direct_exit_initial_placebo_rel_power_tertile2_4",  
  "direct_exit_initial_placebo_rel_power_tertile2_5", 
  
  ## Direct Effect (Power Tertile 3)
  "direct_exit_initial_rel_power_tertile3", 
  "direct_exit_initial_placebo_rel_power_tertile3_1",  
  "direct_exit_initial_placebo_rel_power_tertile3_2",  
  "direct_exit_initial_placebo_rel_power_tertile3_3",  
  "direct_exit_initial_placebo_rel_power_tertile3_4",  
  "direct_exit_initial_placebo_rel_power_tertile3_5", 
  
  ## Indirect Effect
  "indirect_exit_initial", 
  "indirect_exit_initial_placebo_1",  
  "indirect_exit_initial_placebo_2", 
  "indirect_exit_initial_placebo_3", 
  "indirect_exit_initial_placebo_4", 
  "indirect_exit_initial_placebo_5", 
  
  ## Indirect Effect (Absolute Power == 1) 
  "indirect_exit_initial_power", 
  "indirect_exit_initial_placebo_power_1",  
  "indirect_exit_initial_placebo_power_2", 
  "indirect_exit_initial_placebo_power_3",  
  "indirect_exit_initial_placebo_power_4", 
  "indirect_exit_initial_placebo_power_5", 
  
  ## Indirect Effect (Absolute Power == 0) 
  "indirect_exit_initial_not_power", 
  "indirect_exit_initial_placebo_not_power_1",  
  "indirect_exit_initial_placebo_not_power_2", 
  "indirect_exit_initial_placebo_not_power_3",  
  "indirect_exit_initial_placebo_not_power_4", 
  "indirect_exit_initial_placebo_not_power_5", 
  
  ## Indirect Effect (Power Tertile 1) 
  "indirect_exit_initial_power_tertile1", 
  "indirect_exit_initial_placebo_power_tertile1_1", 
  "indirect_exit_initial_placebo_power_tertile1_2", 
  "indirect_exit_initial_placebo_power_tertile1_3", 
  "indirect_exit_initial_placebo_power_tertile1_4",
  "indirect_exit_initial_placebo_power_tertile1_5", 
  
  ## Indirect Effect (Power Tertile 2) 
  "indirect_exit_initial_power_tertile2", 
  "indirect_exit_initial_placebo_power_tertile2_1",  
  "indirect_exit_initial_placebo_power_tertile2_2",  
  "indirect_exit_initial_placebo_power_tertile2_3",  
  "indirect_exit_initial_placebo_power_tertile2_4",  
  "indirect_exit_initial_placebo_power_tertile2_5", 
  
  ## Indirect Effect (Power Tertile 3) 
  "indirect_exit_initial_power_tertile3", 
  "indirect_exit_initial_placebo_power_tertile3_1",  
  "indirect_exit_initial_placebo_power_tertile3_2",  
  "indirect_exit_initial_placebo_power_tertile3_3",  
  "indirect_exit_initial_placebo_power_tertile3_4",  
  "indirect_exit_initial_placebo_power_tertile3_5", 
  
  # Overlap Variables
  "overlap", 
  "overlap_power",
  "overlap_not_power",
  "overlap_rel_power", 
  "overlap_not_rel_power",
  "overlap_power_tertile1", 
  "overlap_power_tertile2", 
  "overlap_power_tertile3", 
  "overlap_rel_power_tertile1", 
  "overlap_rel_power_tertile2",
  "overlap_rel_power_tertile3")

# Rename analysis_dat and save
analysis_dat <- as.data.frame(reform_data[,analysis_dat_vars])
objects_to_save <- c("unilateral_exits", "analysis_dat", "data_full")
save(list = objects_to_save, file = "data/analysis_data.RData")

# Introduction Statistics ----
nrow(regimes) # Total number multilateral treaties
length(unique(analysis_dat$regime_number)) # Number regimes in panel
summary(analysis_dat$year) # Years in time series
sum(analysis_dat$direct_exit_initial) # Number treated regime-years
# See "Substantive Results" section in 03_presentation.R for remaining statistics in introduction

# Data and Measurement Section Statistics ----
length(unique(analysis_dat$regime_number)) # Number regimes in panel
summary(analysis_dat$year) # Years in time series
nrow(regimes) # Total number multilateral treaties
sum(regimes$amendment_code == 1 | regimes$core_treaty == 1) # Total reforms and parent agreements
sum(regimes$amendment_code < 0 & regimes$core_treaty == 0) # Total number of treaties that are not reforms (reported in Appendix A)
sum(regimes$regime_name == "UN Motor Vehicle Regulations", na.rm = T) # Total UN Motor Vehicle treaties
sum(regimes$regime_name == "UN Motor Vehicle Regulations", na.rm = T)/nrow(regimes) # Percent
sum(regimes$amendment_code[regimes$regime_name == "UN Motor Vehicle Regulations"] < 0, na.rm = T) # Total UN Motor Vehicle treaties not reforms
1/mean(analysis_dat$reform_n, na.rm = T) # How often regimes are reformed on average
summary(analysis_dat$reform_n) # Distribution of reforms
nrow(unilateral_exits) # Number of initial unilateral exits
length(unique(data_full$regime_number[data_full$decision_year <= 2018 & # Total number of regimes with unilateral exit
                                    data_full$decision_year >= 1945 & 
                                    data_full$treaty_id %in% unilateral_exits$treaty_id]))
sum(analysis_dat$direct_exit_initial) # Total number of directly treated regime-years
